
> ## ---- Ordered logistic ----
> # Note: postcodes have been pseudonomyzed
> survey <- readRDS("dat/auxiliary/survey.rds")

> ## Estimation
> m <- MASS::polr(as.factor(rent_tri) ~ rent_diff,
+                 data = survey,
+                 weights = weights,
+                 Hess = TRUE)

> vcov_m <- sandwich::vcovCL(m, survey$postcode)

> x <- seq(min(survey$rent_diff), max(survey$rent_diff), length.out = 101L)

> # Parameter simulation
> set.seed(20240513L)

> coef_sim <- MASS::mvrnorm(10000L, c(m$coefficients, m$zeta), vcov_m)

> b <- as.vector(coef_sim[, 1])

> eta <- b %*% t(x)

> k1 <- coef_sim[, 2]

> k2 <- coef_sim[, 3]

> stacked_prob <- abind::abind(t(apply(
+   apply(eta, 2, function (x)
+     plogis(k1 - x)), 2, quantile, c(.5, .025, .975)
+ )),
+ t(apply(
+   apply(eta, 2, function (x)
+     plogis(k2 - x)), 2, quantile, c(.5, .025, .975)
+ )),
+ matrix(1.0, 101L, 3L),
+ along = 3)

> ## Visualization
> pdf(file = "fig/rent_perceptions.pdf",
+     width = 1600 / 288,
+     height = 1200 / 288)

> stacked_prob_plot(
+   est = stacked_prob,
+   x_vals = x,
+   x_lab = "Change in rents, 2015-2020",
+   labels = rev(
+     c("Increased\nstrongly",
+       "Increased",
+       "No increase/\nDecrease")
+   ),
+   offset = -1.5,
+   rug = survey$rent_diff,
+   zeroline = TRUE
+ )

> dev.off()
RStudioGD 
        2 
